function [probabilidade] = prob(p,n,k,N)
%probabilidade e o return output
  % p -> probabilidade de sair cara
  % n -> numero de lancamentos
  % k -> numero de caras pretendidas
  % N -> numero de experiencias
  % PELO MENOS um  (ou seja k == 1 ou k == 2)  -> >=k  (k=1)
  % P(A|B) = P(AB)/P(B)

    lancamentos = (rand(n,N) > p);
    sucessos= sum(lancamentos)==k;
    probabilidade= sum(sucessos)/N;

end

N = 1e5;  %Numero de experiencias
p = 0.5;  %Probabilidade do acontecimento (sair caras)
k = 6;    %Numero de vezes que queremos que saiam caras
n = 15;    %Numero total de lancamentos

probEx2=prob(p,n,k,N)
for k = 0: n
  fprintf("P(Sair %d caras em %d lançamentos) = %d \n",k,n,prob(p,n,k,N));
end


%%analiticamente
prob = factorial(n)/(factorial(n-k)*factorial(k))*p^k*(1-p)^(n-k)

%o histograma representativo da distribuicao de probabilidades
figure(1);
bar(1:6,prob, 'b')

figure(2);
hist(prob,5)